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Optimal control of molecular dynamics is commonly expressed from a quantum mechanical per- 
spective. However, in most contexts the preponderance of molecular dynamics studies utilize clas- 
sical mechanical models. This paper treats laser-driven optimal control of molecular dynamics in a 
classical framework. We consider the objective of steering a molecular system from an initial point 
in phase space to a target point, subject to the dynamic constraint of Hamilton's equations. The 
classical control landscape corresponding to this objective is a functional of the control field, and the 
topology of the landscape is analyzed through its gradient and Hessian with respect to the control. 
Under specific assumptions on the regularity of the control fields, the classical control landscape is 
found to be free of traps that could hinder reaching the objective. The Hessian associated with an 
optimal control field is shown to have finite rank, indicating the presence of an inherent degree of 
robustness to control noise. Extensive numerical simulations are performed, confirming the absence 
of traps in the classical control landscape. Specific cases of controlled anharmonic oscillators are 
presented to illustrate the theoretical principles. We compare the classical formulation with the 
mathematically analogous quantum state-to-state transition probability control landscape. 
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I. INTRODUCTION 



The control of molecular dynamics phenomena with tailored laser pulses is an active research field, and 
the common perspective is to view the dynamics in a quantum mechanical context [T]-[5]. That situation 
is generally necessary when electronic degrees of freedom are involved. However, molecular rotational and 
vibrational motion can often be reliably described classically, as frequently done for treating large polyatomic 
molecules. Some prior works considered optimally controlled molecular dynamics using classical mechanics, 
and in some cases the results suggest that the classical and quantum pictures should have qualitatively similar 
behavior [3141 Oj . The conceptual and computational advantages of classically modeled control warrant further 
analysis of the foundations of the subject. The present paper considers optimal classical control for steering 
a system from an initial point in phase space to a final target point. This work builds on previous classical 
molecular optimal control studies [51 [21 [TTH52] . 

In some contexts, the potential for chaotic behavior of classical molecular systems presents difficulties for 
optimal control. However, some studies have shown that under control, even systems that are nominally 
chaotic may exhibit regular, non-chaotic behavior OUSIHT]. Potential chaotic behavior is most relevant 
in determining the stability of optimally controlled trajectories. In this work, we consider finding such an 
optimal trajectory, but consider its stability only with post-facto bounds on the trajectory encountered. 
Other works, e.g. |^5H26| . have considered using feedback control to stabilize chaotic systems. With chaotic 
systems under open-loop control, the cost function may be re-defined, for example, as an average over a 
window centered around a target time or an ensemble average over initial conditions. An exploration of 
chaotic controlled dynamics deserves special attention, which is beyond the scope of the present work. 

An important feature of optimal control analysis is consideration of the underlying control landscape, 
defined as the physical objective as a functional of the control field. At the critical points of the landscape, 
the functional derivative of the objective with respect to the control is zero. It is important to assess 
whether any critical points correspond to suboptimal traps on the landscape, which might prevent a search 
from reaching the absolute maximum value of the landscape. It has recently been shown that under specified 
conditions the control landscape for maximizing the quantum state-to-state transition probability contains 
no traps [?7^E5] , and the same conclusion holds for more general physical observables [301 131] ■ These latter 
analyses were stimulated by the observed ready achievement of effective control in broad classes of control 
simulations and experiments. In [32], the landscape for optimizing the dynamical transformations in a 
classical system with harmonic potentials is discussed on the symplectic group, and was shown to be devoid 
of traps. A recent paper |33j examined the optimal control landscape of open classical and quantum systems, 
showing that traps do not exist under well-defined conditions. The latter work relies on the open nature of 
the systems, and the present paper considers the dynamics of closed systems. 

In the case of quantum mechanics, prior analyses considered the control landscape for the probability of 
making a transition from one pure quantum state to another — >■ |/)) in a system of N energy levels. 
Here we treat the classical mechanical analogue of making a "pure phase-space state transition" z — ^ z' for 
a system of n particles (i.e., z is a point in the coordinate- momentum phase space). Drawing an analogy 
between a quantum system of N states and a classical system of n particles is unusual, but we show that they 
share strikingly similar landscape topologies. We thus extend previous works that formulate classical and 
quantum mechanics as analogues, either in terms of a classical Hamiltonian system |34[ 135) or Poisson bracket 
operators [35]. The formulations discussed in these papers do not consider the control landscape topology. 
The goal of our work is to (i) analyze the structure of the classical phase space control landscape and (ii) 
perform numerical simulations supporting the theoretical conclusions about the landscape's properties. 



Section jllj introduces the classical system and the target functional J. Section II A considers the control 
landscape critical points where the first-order functional derivative of J with respect to the control field is 
zero. Section [II B| then analyzes the Hessian of the control landscape at the critical points, while Section [II C| 
discusses the robustness and multiplicity of optimal solutions by examining the spectrum of the Hessian. The 
findings in Sections HA and II B are compared in Section HI with analogous quantum mechanical results 
for maximizing pure state transition probabilities. Section ]lV| presents numerical simulations for a simple 
model system to illustrate the results in Section[lT] We discuss simulation results for a one-dimensional Morse 
oscillator and then consider two coupled Morse oscillators. Concluding remarks are given in Section [V] 
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II. OPTIMAL PHASE SPACE CONTROL WITHIN CLASSICAL MECHANICS 



Consider a classically described molecule of n atoms, driven by a linearly polarized electric control field 
£{t). The Hamiltonian H{p, q, t) of the polyatomic molecule can be written as 



^^(P, q, t) = ^P^M-ip + Viq) - D{ci)e{t), 



(1) 



where p = p(t) and q = q(f) are, respectively, 3n-dimensional momentum and position vectors describing 
the state of the n atoms in a three-dimensional Cartesian coordinate frame and 



M 



mi 
m.2 





is a 3n X 3n block diagonal matrix composed of 3 x 3 sub-blocks 



(2) 



m, = 



m, 
m, 

TO, 



; j = l,2,...,n. 



(3) 



with rrii being the mass of the i-th atom. The functions V{q) and D{q) are, respectively, the potential energy 
between the atoms and the projection of the dipole moment along the control field e{t). The analysis below 
assumes that ^(q) and D{q) are twice differentiable functions of the coordinates q and that e{t) is square- 
integrable on the time interval [0,r]. For an arbitrary polyatomic molecule, the 3n Cartesian coordinates 
q may be further written in terms of three center-of-mass coordinates, three overall rotational angular 
coordinates, and 3n — 6 internal coordinates. For simplicity, we do not make the latter transformations, and 
for ease of notation, we set M — 3n in the following treatment. 

The column vector z{t) = [q^{t) P"^(t)] specifies the state of the molecule at time t in the corresponding 
2A/'-dimensional phase space. The dynamical evolution of the molecule is then governed by Hamilton's 
equations: 



dH{p,q,t 

dp 

OH(p.q,t 



/ 

-/ 



Sq 



ag(p.q.t) 

dp 



(4) 



where the dot denotes the time derivative and / is the J\f x J\f identity matrix. Equation Q possesses a 

unique solution [37], given an initial state zq — [qo^ Po"'"]"^- We assume that Eq. ^ has a well-defined 
solution for t S [0,T]. We further make the physically reasonable assumption that the potential and dipole 
functions are bounded, at least over the accessible dynamics on the interval [0,T]. Unbounded potentials 
such as the Coulomb potential can also satisfy our assumptions (e.g., in [38] the Coulomb potential is used to 
classically model hydrogen atom ionization). We assume that the finite time interval implies that all states 
z,(t) driven by admissible control fields are uniformly bounded over t € [0,T] (i.e., the system exhibits no 
finite-time blowups). Thus, an unbounded potential would satisfy our assumption, as long as the points at 
which it is unbounded do not lie within the sampled phase space. The objective is to seek an optimal control 
field to steer the molecule from its initial state Zq to a final state z(T), in order to maximize a specified 
function 0{z{T)). This function is constructed to increase as z(T) comes close to the desired target state 
ztar in phase space. Thus, we seek to assess the maximization of the objective functional J 



maxJ = 0(z(r)) 

e(.) 



(5) 



with respect to the control field e{t), t e [0,T]. The control landscape J = J[e(t)] is a functional of the 
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control field. We assume that all final states z(T), including those that maximize 0(z(T)), lie in the interior 
of the reachable set from Zq at time T. The latter set consists of all final states at time T to which the initial 
state Zq may be steered for some admissible control field e{t) [27-30] /~ Moreover, for simplicity, 0(z(T)) is 
assu med to be twice continuously differentiable with respect to z(T) and bounded from above (see Section 



IIBl. Other terms may be included in the objective functional J (e.g., J = 0{z{T))+a e{ty dt, a > 0), to 
bias the class of controls. As in the quantum mechanical analysis [SHlES], here we focus on the fundamental 



case where J[e(t)] = 0(z(T)) to assess the landscape topology without additional objectives or constraints 
on the field. Choosing 0{z{T)) so that it rises as z(T) — >• Ztar does not in itself assure that J is trap-free 
with respect to e(t). Analyzing the behavior of J with respect to e(t) is an important goal of this paper. 



A. Critical points of the control landscape 

The critical points of the control landscape J[e{t)] correspond to the zero gradient condition 



SJ 

W) 



dJ 6z.{T) 
dz{T) 5e{t) 



= o,yte [0,T]. 



(6) 



— Sz{T) is of full rank (often 
A regular control e(t) is thus 
(5e(t) can move the final state 

the vector-valued function of 



The control landscape critical points may be divided into two types: kinematic critical points, at which 
dJ/dz{T) — 0, and non-kinematic critical points, at which dJ/dz{T) 7^ 42 . Kinematic critical points 
may be further considered as regular, at which the functional mapping de{-) 
referred to as being surjective), or singular, at which is not of full rank |37j. 
defined as having the property that an appropriate perturbation e{t) e{t) -\ 
z(T) — >■ z(T) -I- Sz(T) in any specified direction. 

We assume that the functional mapping 6e{-) — > 5z{T) is of full rank (i.e. 
time Sz{T) / 6e{t) is of full rank 2Af (= 6n)) for the field e{t) at the critical point of the control landscape 
under consideration. Using Eq. ^ it then follows that dJ/dz{T) = at a critical point. Numerical 
simulations (see Section IV I suggest that singular controls are rare, which is consistent with the assumption 
of typical fields being regular and 5z{T)/5e(t) being surjective. We will also assess the likely existence of 
regular controls through an explicit expression derived below for Sz{t)/Se{t). A bound on 6z{t)/Se{t) will 
also be identified, which will be used in Section |II C| to explore the robustness of the control solutions to 
noise at the critical points. 

Consider a perturbation of the control e(t) — ^ e(0 + ^^i^) Eq. Q [39] with the corresponding response 
in z{t) — z{t) + 5z{t) leading to the first order variation 



5z{t) = A{t)6z{t) + B{t)6e(t), 



(7) 



where Sz^{t) ^ [Sq^it) dp^it)]' 



Ait) 



)V(q) 
0q2 



aq2 



M-1 



(8) 



and 



Bit) = 



aj'(q) 

Oq 



(9) 



Here is the J\f x J\f square zero matrix in A{t), and the first TV entries of B{t) are zero. The notation 
d^X/dq^ in A{t) refers to the following J\f x J\f second derivative matrix of the function X{q) (i.e., for ^(q) 
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or D(q)): 



d_ /ax 

aq \ 5q 



dqf 
d^X 



d^X 



d-'X 
dqidq2 



d-'X 



dqj^dqi dqMdq2 



d-X ' 
dqidqM 

dqidqj^ 



d^X 
9ll 



(10) 



Similarly, i9D(q)/9q is the A/'-dimensional gradient row vector of -D(q) with respect to the N components 
of q. Integrating Eq. Q yields 



^z(r) = M[T) / Nr^{t)B{t)5e{t) dt, 
Jo 



where M{t) satisfies the equation 



M{t) = A(t)M(t), M{0) = /. 



(11) 



(12) 



The elements of M have the meaning Mij{t) = dzi{t)/dzj{0), representing the sensitivity of the ith compo- 
nent of the state at time t to the jth component of the state at the initial time. The matrix M{t) is also 
symplectic [44 and thus invertible. From Eq. (11) we have 



Se{t) 



= M{T)M-'^{t)B{t). 



(13) 



Based on Eq. (13), the regularity criterion for a field e(t) or equivalently the surjectivity of 6z{T)/Se{t) 
corresponds to the condition that the vectors M^^{t)B{t) span M.'^-^ , since M{T) is an invertible matrix. 
In the Appendix, we explicitly demonstrate regularity for critical point controls of a one-dimensional linear 
forced harmonic oscillator, for which the vector M~^{t)B{t) is of full rank (i.e., of rank 2) even though the 
corresponding matrix A{t) in Eq. ( |12[ ) is independent of time t. For a general n-particle Hamiltonian system, 
the matrix A{t) in Eq. ([s]) will likely be a complicated function of time (via the state variables q{t) and 
the field e{t)). It appears difficult to prove surjectivity under all conditions, but for arbitrary Hamiltonians 
and fields the mere complexity of the dynamics suggest that the vector components of M~-^{t)B{t), as a 
function of time < i < T, are likely of full rank 2Af. We adopt this assumption here and investigate its 
consequences, which may be tested in numerical simulations. 



In Section [nBj we use the second-order functional derivative of J with respect to the control to investigate 
the optimality of regular critical points and the associated robustness to field errors e{t) — >■ e{t) + 5€{t). 
These investigations enable a post-facto evaluation of the stability of the evolved system trajectory. The 
latter analysis utilizes the time-averaged norm of the functional derivative Sz{T)/Se{t): 



where 



fe(T) 



Se{t) 



\B{t)\\i dt = 



T 



dq 



dD{q) 
dq 



2 

T ' 



dt. 



(14) 



(15) 



which is expected to be finite for typical molecular dipole moments and bounded dynamics q(t), < < < T. 
Here || • II2 is the Euclidean norm of a vector and |j • Ijp is the Frobenius norm of a square matrix. Moreover, 
we have assumed that \\M{T)M-'^{t)\\^ is uniformly bounded (i.e., \\M{T)M-^{t)\\^ < C for aU t £ [0,T], 
C being a finite positive number). The entries of M^^{t) may be bounded by exponentiating a uniform 



bound on the entries of A{t) using Eq. (12). 
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From Eq. the existence of such a bound on the entries of A corresponds to a uniform bound on the 
second derivatives of the potential and dipole moment functions. Since we assume that the state trajectories 
are bounded and that ^(q) and D{q) are twice continuously differentiable, an upper bound on the entries 
of M~^{t) exists. Furthermore, we can then show that the time average of the gradient SJ/de{t) is bounded 
as follows: 



SJ 



Se(-) 



dJ (5z(T) 



9z(r) 5e{t) 



dt < 



C 



dJ 



dziT) 



\B{-)\\, 



(16) 



The result in Eq. 



( 16 ) shows that extraordinarily steep regions should not exist on the landscape, thereby 

Given these critical point conditions and 



favoring the stability of algorithms seeking to optimize J[e(t)] 
gradient bounds, below we consider the optimality of regular critical point controls, using the second-order 
functional derivative of J with respect to the control field. 



B. Nature of the control landscape critical points: Hessian analysis 



This section considers the second-order functional derivative H(i,<') = S^J/Se{t)Se{t') (i.e., the Hessian of 
the objective with respect to the control) evaluated at the critical points. The properties of Hit, t') determine 
the conditions under which the landscape may contain traps (i.e., controls producing local, suboptimal 
maxima). Using the relation 6J/6e{t) — [dJ/dz{T)] [Sz{T)/Se{t)] and differentiating again, the Hessian 
T-L{t,t') can be written as 

^' ' dz{T) 5e{t)5e{t') \5e{t') ) dz{T)^ 6e{t) ' ^ ' 

where 

d^J _ d f dJ 



dz{TY ~ dz{T) \dz{T)^ 
At kinematic critical points, dJ/dz(T) — (see Section II A[ ) and the Hessian becomes 

y(t t') = ^'-^ = f Sz{T)Y J^SziTl 

^ ' ' Se{t)Se{t') \6e{t') J dz{T)^ Se{t) ' ^ ' 

which is symmetric and separable, thus at most of finite rank 2Af, dictated by d^J/dz{T)^. Here we use the 
assumption that J is twice continuously differentiable with respect to z{T) to conclude that d'^J/dz{T)^ is 
a symmetric matrix. 

Establishing the positive, negative or indefinite nature of the Hessian is essential for assessing whether 
traps exist on the landscape. We can test for the character of the Hessian by diagonalizing the matrix 
J / dz{T)^ ; since this matrix is real and symmetric, there exists an orthogonal matrix P such that 

8"^ J 

P^SP, (19) 



dz{T) 

where S is a diagonal matrix, whose diagonal entries tXi, i = 1,2,..., 2A/' are the eigenvalues of J jdziT)^ 
[15] . Then Eq. (|l8]) becomes 
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which may be re-written as 



2A/' 



(21) 



1=1 



where each (j)i{t) is the Ith component of the vector P {Sz{T)/Se{t)). Thus, the rank of the Hessian H{t,t') 
equals the number of nonzero eigenvalues (Ti of J / dz{T)^ . Moreover, it can be shown [IB] that the signature 
of the Hessian (the difference between the number of pos itive and number of negative eigenvalues) is the 
signature of S, or equivalently 9^J/9z(r)^, from Eq. (19). 

Since J = 0(z(T)), we see that a regular control field e(t) corresponds to a local or global maximum if 
it produces a final state z{T) such that d^J/dz{T)'^ is negative definite. Importantly, we are free to choose 
the observable J — 0(z(T)) on physical grounds so that it contains only one maximum which is global, and 
in this case the control landscape contains no traps for regular control fields. Moreover, the landscape will 
contain saddle points at regular control fields if and only if 0(z(T)) has saddle points. Thus, traps on the 
landscape under the assumptions of the analysis should only arise if they were artificially introduced into 
0(z(T)). 

Establishing a bound on the Hessian is important (see Section al C ), as it reflects the degree of robustness 
at the optimal value of J to variations in the field e{t) — )■ e{t) + 6e{t) due to noise. We can obtain a bound 
on the trace of the Hessian 'H{t,t') using Ec^. (pO|: 



iTr-Hl 



'H{t,t) dt 



< C^T\\B{- 



1M 
1=1 



(22) 



For example, if 0(z(T)) is the quadratic function [z(r) — ztar]"^W[z(r) — ztar], where W is a symmetric, 
negative-definite 2A/'x IM matrix, then 0{z{T)) has only one critical value (the global maximum) 0{z{T)) = 
0, corresponding to z{T) — ztar = 0. In this special case, it is readily seen that the kinematic Hessian at the 
global maximum is the 2J\f x 2J\f matrix J / dz{T)^ = 2W, which is negative definite and of full rank 2Af 
by assumption. Therefore TrV. will be negative, as Tr {d'^J/dz{Tf) = 2Tr W < 0. 



C. Robustness analysis and level sets at the control landscape global maximum 



The existence of only a finite number of nonzero eigenvalues for 'H{t,t') at controls producing global 
maxima implies that the Hessian has infinitely many zero eigenvalues (i.e., T-Lit, t') has an infinite-dimensional 
nullspace). This property has important consequences for the robustness to noise and the multiplicity of 
optimal controls. Qualitatively, a control solution e{t) is robust if a perturbation 5e{t) does not significantly 
change the value of J . The fact that the Hessian possesses infinitely many zero eigenvalues indicates that 
the landscape is inherently fiat at the global maximum value with the existence of a submanifold of optimal 
control solutions around a given optimal field. 

The robustness of an optimal control solution follows from the previous bounds on the trace of the Hessian 
in Eq. (22). Consider the response SJ due to a perturbation 5e{t) up to second-order 



5J 



6J 

W) 



8e{t) dt 







5'J 



5e{t)Ht') 



5e{t)6e{t') dt' dt. 



Using Eq. (21 1 at the control maximum (i.e., SJ/Se{t) — 0), Eq. (23) implies that 

, 27V / T \^ 7^2 27V 

M = 2^'^'' / / Mt)Mt')Se{t)Se{t')dtdt'\ < -WSef^ x Y,W,\\\M 

1 = 1 -'^ J 1 = 1 



2 



(23) 



(24) 



where we used u/ < 0; ^ = 1, 2, . . . , 2Af for a control corresponding to the landscape's maximum value. The 
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norm |l(5e||T is defined as in Eq. (15): 



i / Mt)f dt. 



(25) 



Since the trace of the Hessian Tr H satisfies 



•^0 ;=i -^0 1=1 



Eq. (24) can be rewritten as 



\5J\ < ^WSefr X |TrH| < Jc^T^H^elH x \\B{-)\\1 x 



Tr 



_(PJ_ 



(26) 



(27) 



where we used the fact that Tr (^d^ J / dz{T)^^ < 0. This result shows that a near-optimal control field (i.e., 
e{t) + Se{t)) will produce a near-optimal objective functional value J. The bound on the resultant variation 
SJ scales with the squared- norm of the control field and the magnitude of the trace of the second derivative 
of J with respect to the final state. This result along with the infinite-dimensional nullspace of H{t,t') 
implies the existence of inherent robustness to control noise at the landscape maximum value. 



The analysis above also suggests that a level set of optimal controls exists near any given optimal solution 
e(i), since the objective functional J remains maximized when the control e is systematically altered by a 
variation that lies in the infinite-dimensional nullspace ol'H{t,t'). To move the control in this nullspace, we 
adapt the procedure exploited for the analogous behavior of quantum systems [47 . The control field is now 
denoted as e{s,t), where the variable s > labels a point in the evolution of the control field within the 
nullspace of H{t,t'). Since the gradient satisfies SJ/6e{s,t) = at the optimum, the first derivative dJ/ds 
is automatically zero. 



dJ_ 
ds 



5e{t) ds - ' 



and we wish to choose de{s,t)/ds so that the second derivative satisfies d^J/ds^ = 0: 



^ '''^5e(s,0^(^^^„ae(s,i) 



ds 



ds 



2M 

dt dt' — ai 



de{s,t) 



-\ 2 



dt 



= 0. 



(28) 



Thus, we require that for some arbitrary, i^-integrable function f{s,t), 



de{s,t) 
ds 



2Ar „T 

^f{s,t)-y2Mt) / Mnfis,ndt', s>o 

1 = 1 -^o 



(29) 



at the global maximum, where the functions ^i{t) are orthonormal (with respect to the || • \\t norm) with 
the same span as the functions 4>i{t) in Eq. ( [2l] ). The functions ipiii) may be obtained by applying the 
Gram-Schmidt process to the functions 0i(i). Equation (29) automatically satisfies d'^Jjds^ = 0, which is 
evident upon direct substitution into Eq. 



( p8| . The differential equations in (29) are coupled to Hamilton's 
equations (|4]) and initiated with an optimal field e(0, i) already found to produce the global maximum value 
for J. 
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III. COMPARING THE CLASSICAL AND QUANTUM CONTROL LANDSCAPES 



Our analysis of the classical phase space trajectory control landscape has a mathematical analogue with the 
quantum state-to-state transition probability control landscape [28, i29i . In particular, optimally controlling 
the 27V-dimensional quantum wavefunction (i.e., its joint N real and N imaginary components) of a system 
with TV energy levels is mathematically similar to optimally controlling the 2Af position and momentum state 
variables of n classical particles. We emphasize that this analogy is not a direct physical comparison of the 
two classes of systems. The latter distinction is evident, as controlling the motion of an n-particle quantum 
mechanical system requires a wavefunction which lies in an infinite-dimensional space. Nevertheless, there are 
striking mathematical similarities between the n-particle classical and iV-state quantum control landscapes, 
as these two cases describe analogous state-to-state "transitions" in their respective classical and quantum 
regimes. 

The classical analogue of the quantum transition probability is the corresponding observable function 
0(z(T)) = [z(r) — Ztar]"^W[z(T) — Ztar], witli W Symmetric and negative-definite and ztar the target state. 
Since 0{z{T)) has only one critical point, at which z(T) = Ztar and 0(z(T)) is globally maximized, the 
analysis in Section |ll] shows that the classical target-state landscape is trap-free with no saddle points. 
To obtain this result, we assumed that the functional mapping i5e(-) — )■ Sz{T) is surjective, implying that 
appropriate perturbations Se{t) in the control field can move the final state z(T) in an arbitrary direction 
6z(T). Quantum mechanical studies make the analogous assumption that proper perturbations in the control 
field can enable the final state \ip{T)) to vary in any direction |(5'i/;(T)) [251 [^ . 

In order to show the analogue of the quantum mechanical case to the classical formulation in Section [ll] we 
summarize here the key quantum mechanical state-to-state landscape relations. The objective functional for 
quantum mechanical state-to-state transition probability control may be written as Jqu = {■ip{T)\f) (/|V'(T)), 
where |V'(i)) evolves from the initial state \i) at i = to the final state {ipiT)) [351[5H]. The goal is to maximize 
Jqu (i.e., steer \ip{T)) to be aligned with |/), up to an overall phase) in analogy with the classical goal of 
maximizing J in Eq. ([s]). As with the critical point condition for the classical control landscape (Eq. ([6])), 
the critical point of the quantum control landscape Jqu[e] corresponds to the condition that the gradient 
5Jqu/Se{t) be zero: 

Seit) UlV'(r))' Seit) 

where we used the convenient notation {w, v) — Re (w^v'j to express the inner product of two iV-dimensional 
complex vectors. Here "Re" is the real part of a complex number. 

The evolution of the wave function \4>{t)) follows the time-dependent Schrodinger equation {d/dt)\i{}{t)) — 
{l/{ih))(^Hf) — ne{t)^\il!{t)), \ip{0)) = \i), where Hq and /i are, respectively, the field-free Hamiltonian and the 
dipole moment operator of the quantum system. In order to avoid confusion with the classical Hamiltonian 
m Eq. ^, we use explicitly different notation for the latter quantum mechanical operators. For an TV-level 
quantum system, due to (i) the preservation of the norm {ijj{t)\ifj{t)) = 1 and (ii) the fact that the overall 
phase of \ip{T)) is not physically relevant, the state \tp(t)) may be considered as a column vector of up to 
2N ~ 2 independent components. 

The propagator f7(i,0) satisfies the equation {d/dt)U{t,0) = -{i/h) {Hq - fJ,e{t)) U{t,0) with J7(0,0) = 
I, which is the quantum analogue of M{t) in Eq. ( |l2| ). We can solve for the functional derivative 
5\tlj{T))/6€{t) = {i/h)U{T,0)W{t,0)fiU{t,0)\iP{0)) and use the relation dJqJd\^/j{T)) = 2|/) (/|V'(T)) to 
obtain the gradient 

^ = (2C/t(T,0)|/)(/|^(r)),^t/t(t,0)M{7(t, 0)1^(0))), (31) 

which leads to critical point controls when Jqu = or Jqu — I, respectively, at the global minimum 
((/IV'(T)) = 0) or at the global maximum (J7l'(r, 0)|/)(/|f/(r, 0) = \i){i\). At both of these "regular" 
critical points, S\iIj{T)) /Se{t) is of maximum 27V - 2 rank [2511^ . 
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At the global maximum the corresponding Hessian can be written as |27) 



'Hqu(i, t') 



Re 



(32) 



which is negative-semidefinite and at most of rank 2N — 2. The summation in Eq. (32) is analogous to that 
of Eq. pTj ); both are sums of linearly independent functions, with the number of such functions equal to the 
rank of the derivative of the state with respect to the control. Here we may associate \/\cri\4>i{t) in Eq. (21 ) 
with {k\ii{t)\i) in Eq. (32). Thus, as with the classical circumstance, the Hessian for the quantum state- 



to-state transition probability is symmetric and of finite rank when we assume surjectivity of the functional 
derivative S\'ip{T)) /5e{t). In addition, the norm of 6Jqu/S£{t) and the trace of the quantum Hessian are 
bounded [SSj, just as for their classical counterparts (see Section HC). Moreover, one may explore the level 
sets of quantum landscapes, much as Section |H C| considered the level set at the global maximum in the 
classical case [iSl. 



IV. NUMERICAL SIMULATIONS 



The results in Section|lT]on the existence of trap-free classical state-to-state control landscapes rest on some 
key assumptions, particularly the surjectivity of 6z{T)/Se{t). To assess to validity of these assumptions, we 
performed classical optimal control calculations based on Hamilton's equations for a range of model molecular 
systems with various Hamiltonians. In these calculations, we considered the objective functional 



J = - [(q(T) - qtar) (p(r) - ptar)] ^ K 



(q(T) - qtar) 
(p(T) - ptar) 



- -fcl(q(T) - qtar)^(q(r) - qtar) - fc2(p(T) - Ptar)^(p(T) - Ptar), (33) 

where (q(r), p(T)) and (qtar, Ptar) are pairs of real vectors denoting, respectively, the final and target states 
of the molecule in phase space. The matrix K = diag (fci, ^2), with fci > and k2 > scaling the position 
and momentum targets; unless otherwise specified, we take K to be the identity matrix (fci =^2 = 1)- 
The optimal control simulations were carried out using fourth-order symplectic integration of Hamilton's 
equations [35] in conjunction with a gradient-based algorithm analogous to the D-MORPH method, which 
has been utilized in quantum optimal control simulations |47[ I50| . Many simulations were run, and a few 
cases will be shown to illustrate the principles of Section [H] We first present two cases for simple vibrational 
control of a heteronuclear diatomic molecule aligned along the field direction. We then present the results 
of several simulations with systems of coupled oscillators and 2n = 4 position and momentum coordinates. 
In all of the tests, no traps on the landscape were encountered. 



A. Numerical techniques: Symplectic integration and D-MORPH optimization algorithms 

In the following simulations, the initial trial control field eo(t) is chosen as a sum of four sine functions 
with randomly selected amplitudes and phases. The field is allowed to freely vary as a discretized function 
of time, < i < T, during the optimization. The control field e{t) evolving towards its optimal form is 
discretized on an evenly spaced temporal grid, i.e., tj — jAt e [0, T], j = 1,2,..., [T/ At\ , where At is a 
small, fixed time-step. Finding an optimal control field involves iterating two steps: (i) solving Hamilton's 
equations and (ii) updating the control field, until J has reached a maximum value. 

In step (i), a fourth-order symplectic integrator [5T] was used to integrate Hamilton's equations (|4| for 
the state evolution over the discretized time interval [0,T]. A symplectic integrator was chosen for reasons 
of numerical stability, although other methods could be utilized. The adopted method uses four substeps 
to evolve the state from time t to time t + At. Specifically, denoting the state at time t by (q°,p°), for 
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2 = 1,2, 3, 4, we have 



OH , , 1 . . 



(34) 
(35) 



with {q{t + At),p{t + At)) = (q^'iP^). The vahies of the fixed parameters bi and Ci are given in [491152]. 

In step (ii), the control field is refined using a gradient search algorithm based on D-MORPH [53] via a 
homotopy parameter s, such that the field is written as e{s,t) and e{0,t) — eo(i). The D-MORPH equation 
to update the field is 



de{s, t) 
ds 



/3 



6J 



6e{s,t) 



/3>0, 



ensuring that 



dJ 

ds 



5 J de{s,t) 
Se{s,t) ds 



dt^ /3 



6J 



(5e(s,i) 



dt > 0. 



(36) 



(37) 



Equation (36 1 is discretized over time, as described above. The field e{s,t) evolves as £{s,t) — e(s + 
As,t), thereby increasing J. Integration of Eq. (36) is stopped at the convergence criterion: —J < 0.01. 
Equation ( 36 ) is integrated with respect to s at each of the time points tj using a fourth-order Runge- 
Kutta method. At each increment in s, the gradient SJ/Se{s,t) is computed with Eq. ([6| and dJ/Se{t) 



[dJ/dz{T)] [Sz{T)/Se{s,t)] once Hamilton's equations are solved, using Eqs. (34) and (|35[) in step (i) based 
on the control field in the preceding iteration. 



B. Example: Optimal control of a diatomic molecule 



Consider a model vibrating diatomic molecule driven by a time-dependent control field e(t), t G [0,r], 
aligned with the molecular axis. The diatom is modeled as a Morse oscillator with Hamiltonian 

- y + V{q) - D{q)e{t) ^ ^ + Do {l ~ e'^'^f - Aqe-^'^\{t), (38) 



where q and p = p/^/m are, respectively, the relative distance and scaled momentum between the two atoms. 
The center of mass motion has been removed and the problem is expressed in internal coordinates; however, 
the control landscape in this formulation should still exhibit the general behavior identified in Section [III 
For illustration, the parameters for the Morse potential V{q) and dipole moment function D(q) in Eq. 



( |38[ ) are chosen to have the values: m = 1732, Dq = 0.2101, a = 1.22, A = 0.4541 and ^ = 0.0064, all given 
in atomic units (a.u.) |T5l [54l [55] . We define e{t) — G{t)E{t), where G{t) is a Gaussian envelope function, 
and then optimize with respect to E{t). The presence of the Gaussian guides the fields e{t) towards the 
physically desirable behavior of being small at the beginning and end of the time interval [0, T]. Two cases 
will be shown here, operating respectively in the strong and weak field regimes, to particularly illustrate that 
the principles in Section [H] are expected to be broadly applicable. 

In the first case, the initial condition is chosen as (gojPo) = (O-^: 0-0), corresponding to a slightly stretched 
bond (i.e., the equilibrium position is g = 0) with zero momentum. The desired target state at T = 3207r 
is (gtanPtar) = (2.0,0.014), which is far from the initial state. Note that for the momentum, the distance 
between the initial value of 0.0 and the target of 0.014 is quite substantial, as we will demand precision at this 
level and the phase space excursion will range out to \p\ 0.8. As the time interval is also relatively short, 
we expect that the optimal field will be strong. For numerical stability, we take K — diag (0.999, 1.732) in 
Eq. 

Figure [l] shows the initial and final optimal control fields. The objective had the value J = —3.25 for the 
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FIG. 1. Initial field and final optimal control for the single oscillator with a strong field. The initial field yields the 
phase space point {q{T),p{T)) = (0.19,3.6 x 10~^). The target state was (gtar, Pj^^^,) = (2.0,0.014), and the optimal 
field achieved {q{T),p{T)) = (2.000,0.015). 



-Initial 
-Optimal 




Position (a.u.) 



FIG. 2. Phase plane trajectories with the initial and final optimal control fields for the single oscillator. The point 
denoted with x marks the final state of the trajectory with the optimal field, while that denoted with + marks the 
final state with the initial field. The square ■ marks the initial point {qo,pg). 



initial field, and the final field achieved the near optimal value J = —1.87 x 10 corresponding to the final 
state q{T) = 2.000 and p{T) = 0.015. The latter values are very close to the target state of (2.0, 0.014). 

Figure |2] shows the phase plane trajectories arising from the initial field and final optimal field. The 
optimal control field substantially changes and "stretches out" the q trajectory in order to reach gtar = 2 at 
the final time. This behavior is consistent with the optimal field having a large amplitude towards the end 
of the trajectory, as shown in Fig. [T] 

A key feature of the analysis in Section [ll] is the assumed surjectivity of 6z{T)/6e{t). For the illustration 
here, we calculate the singular values of the time-discretized matrix 6z{T) /de{t) to be 1.77 and 0.01, showing 
that the two functions 6q{T) / 5e{t) and 6p{T)/6e{t) are comfortably linear independent. 

We next show an example of the same diatomic model with a weak optimal field. The initial and target 
states are chosen to be (go,Po) = 0.012) and (gtarjPtar) = (O-^i 0.12), respectively, and we substantially 

lengthen the control interval to T = 32007r. The matrix K in Eq. (|33]) is again taken to be diag (0.999, 1.732). 

Figure [3] shows the optimal control field and resulting phase plane dynamics. We see that the field is much 
weaker than that in Fig. [T] with a maximum amplitude less than 0.06 a.u. The phase plot dynamics are 
quasi-periodic, and come quite close to the target state: with the optimal field, {q{T),p{T)) — (0.497,0.119), 
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(a)Optimal control field. (b)Pliase plane dynamics. 

FIG. 3. Optimal control field and phase plot for the single oscillator with weak control field. The 'x' marks the final 
point with the optimal field, and the marks the final point with the initial field; the square □ marks the initial 
point (go,Po) = (-0.3,0.012). With the initial field, the final state is {q{T),p{T)) = (-0.300,-0.171) and with the 
optimal field it is {q{T),p{T)) = (0.497,0.119). The target state is (gtar,Pta,r) = (0.5,0.120). 



-2 




1 2 3 4 5 6 7 

S 

FIG. 4. Evolution of the objective functional J versus s, the generalized iteration number, for the single oscillator 
with weak control field. The optimization was stopped when J = —8.46 x 10~^. 

for an objective functional value J = —8.46 x 10^^. Figure |4] shows the corresponding evolution of the cost 
functional. 

Finally, we examine the surjectivity of Sz{T) / Se{t) at the optimal field. For this illustration, Fig. [5] shows 
the two functions 5q{T) / 5e{t) and Sp{T) / Se{t) with the optimal field. We see that these are visibly linearly 
independent; the singular values of the time-discretized dz{T) / 6e{t) matrix are 3.978 and 0.236, showing a 
comfortable linear independence. 



C. Example: Optimal control of two coupled Morse oscillators 

We consider a two-dimensional coUinear system of coupled Morse oscillators, as examined in [S5J[57]. We 
choose a coordinate frame where the two oscillators are coupled through the potential and the momenta 
are decoupled. The system is driven by a time-dependent control field e{t), t £ [0,r]. We let {qi,q2) and 
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FIG. 5. Presentation of 5z{T)/Se{t) for the single oscillator with weak control field. The singular values are 3.978 
and 0.236, thereby showing clear linear independence of the two functions. 



(Pi)P2) — (Pi/%/"^iiP2/\/'Ti2) denote the two position and scaled momentum coordinates respectively. The 
Hamiltonian is given by 

H=l (pI+pI) + Viqi,q2) - D{q,,qMt) 



—9 —9 

Pi pi 

'22 



(1 - e-^'i^y + (1 - e-"'"'^y + 2^ (1 - e-"!"!) (l - e""^"?^) 



q2e-^''^)e{t). 



(39) 



The center of mass motion has been removed and the problem is expressed in internal coordinates; as in the 
single oscillator example, the control landscape in this formulation should still exhibit the general behavior 
identified in Section HIl 



The parameters for the potential ^(91,92) and dipole moment function D{qi,q2) in Eq. (39) are chosen 
to have the values: Dq = 0.2102, ai = 1.1282, a2 = 0.9712, /? = 0.1, and ^ = 0.1, all given in atomic units 
(a.u.) 56 . The masses are chosen as mi = 10.0, m2 = 5.0. Sixteen distinct simulations were performed with 
the initial conditions for qi, q2, Pi, P2 chosen randomly at intervals of spacing 0.2 between —1 and 1 for the 
position coordinates; the momentum coordinates Pi and P2 were chosen at intervals of 0.06 between —0.32 
and 0.32 and 0.09 between —0.45 and 0.45 respectively. The equilibrium position of the field-free Hamiltonian 
is at qi ^ q2 = 0. The desired target state at T = 207r is (gtar.ij 9tar,2jPtar nPtar 2) = (2-0, 0.5, 0, —0.22). 

Figure [6] shows the increasing objective (i.e., log(— J)) as e(s,t) evolves with s for each of the sixteen 
simulations; here J is given by Eq. (33 1 with K = diag (1, 1, 0.1, 0.2). The initial control field was chosen as 



a sum of four sine functions with random amplitudes between and 1 and frequencies at approximately 0.5, 
1.0, 1.5, and 2.0 times the frequencies of the (gi,Pi) and ((72, P2) field-free motion. The frequencies of the two 
oscillators are close to each other due to the choices of ai and a2 in the Hamiltonian ( 39 ) . The increase in J 



generally is very rapid at the beginning and slows down as s becomes large, due to the diminishing gradient 
SJ/6e{s, t). In some cases (e.g., simulations j and p), the rate of increase is initially slow, and then increases. 
These occurrences likely correspond to the evolution starting out at areas in the landscape of lower slope. 
Nevertheless, the increase in the objective functional was strictly monotonic in all cases with no evidence of 
traps. No attempt was made to traverse the level set at the maximum of J as discussed in Section [H] The 
final state coordinates at the optimal control field are shown in Table |l] 

Figure [7] shows the initial and final optimal control fields for the simulation o; no attempt was made to 
impose the constraint that the field die out at ^ = and T. The moderate difference between the initial and 
final fields in Fig. [7] has a dramatic impact on the final achieved state as shown in the figure caption. The 
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FIG. 6. The evolution of the objective functional J versus s for sixteen different simulations using the two coupled 
oscillators. The objective functional value approached zero in all cases, and the optimization was stopped when the 
values reached those shown in Table [I] 



Simulation 


91 (0) gi(T) 


92(0) 92(7) 


Pi(o) Pi(r) 


P2(0) MT) 


a 


-0.2 


1 


982 


0.2 


0.504 


-0.19 


0.007 


-0.18 


-0.225 


b 


0.8 


1 


958 


0.8 


0.495 


0.19 


0.013 


-0.27 


-0.213 


c 


0.2 


1 


967 





0.491 


0.13 


0.013 


0.18 


-0.211 


d 


0.4 


1 


957 


1.0 


0.493 





0.013 


-0.18 


-0.210 


e 


-0.2 


1 


956 





0.493 





0.014 


0.36 


-0.209 


f 


0.4 


1 


948 


-1.0 


0.491 


-0.19 


0.016 


-0.27 


-0.207 


g 


1.0 


1 


964 


-1.0 


0.495 





0.011 


-0.09 


-0.212 


h 





1 


962 


-0.8 


0.491 


-0.13 


0.014 


0.09 


-0.209 


i 


-0.8 


2 


018 


0.2 


0.498 


0.13 


-0.008 


0.09 


-0.223 


j 


-0.6 


2 


015 


-0.8 


0.508 


0.06 


-0.007 


-0.09 


-0.230 


k 


-0.6 


2 


025 


0.2 


0.500 


0.06 


-0.012 





-0.225 


1 





1 


957 


0.4 


0.494 





0.014 


0.18 


-0.208 


m 


-0.2 


2 


019 


0.6 


0.511 


0.06 


-0.009 


0.27 


-0.244 


n 


0.8 


1 


976 





0.494 


-0.06 


0.009 


0.18 


-0.215 





-0.6 


2 


023 


-0.2 


0.489 


0.13 


-0.011 


-0.18 


-0.218 


P 


-0.6 


2 


016 


-0.6 


0.503 


0.25 


-0.009 


-0.09 


-0.225 



TABLE I. Initial state conditions and final state coordinates with the optimal fields of the simulations shown in Fig. 
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FIG. 7. Initial and final optimal control fields for simulation case o with the two coupled oscillators. The initial 
field yields the state coordinates {qi{T), q2{T),p^{T),P2{T)) = (10.61, -0.40, 0.52, 0.23). The target state was 
(qtar.i,Qtar.2,Ptar,i,Ptar.2) = (2.0,0.5,0,-0.22), and the Optimal field achieved {qi{T),q2{T),p,iT),p^{T)) = (2.02, 
0.49, -0.013, -0.22). 




(a)Position. (b)Momcntum. 



FIG. 8. Position and momentum trajectories for case o of the two coupled oscillators with initial and final control 
fields shown in Fig. [7] 



objective had the value J — —78.71 for the initial field, and the final field achieved the much higher and near 
optimal value J = -0.0019, corresponding to the final state qi{T) = 2.02, g2(T) = 0.49, Pi(r) = -0.013 
and and P2{T) — —0.22. The latter values are very close to the target state of (2.0, 0.5, 0, -0.22). Figurejs] 
shows the state trajectories for the initial and final optimal control fields in this example trajectory. The first 
position coordinate qi{T) with the initial field has a value far above the target of 2.0. The optimal control 
field substantially changes and "compresses" this qi trajectory in order to reach gi(T) = 2.02 at the final 
time. Similar behavior was found for the other position and momenta coordinates. The assumed surjectivity 
of (5z(T)/i5e(s, i) was also checked, and Fig. [9] shows the eigenvalues of the time-discretized Sz{T)/Se{s,t) 
matrix as the control field e evolves with s. The eigenvalues show comfortable linear independence of the 
rows of dz{T)/Se{s, t) for all s values. Similar results were observed for the other examples shown in Fig. [6] 



V. CONCLUSION 



This paper examines optimal control of a single phase space trajectory for a system following Hamiltonian 
classical dynamics. We have shown that for these systems, all critical points of the optimal control landscape 
J[e(i)] are critical points of the observable 0(z(r)), under specified assumptions including surjectivity of 
the functional mapping 5e{-) Sz(T) and boundedness of the potential function over the sampled phase 
space. In particular, if the observable is a quadratic function aiming to reach a target state, a gradient-based 
search starting from an arbitrary field eo(i) should reach a critical (optimal) field e(t) achieving this target 
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FIG. 9. The singular values of 5z(T)/5e(s, as a function of s over the control field evolution to climb the landscape 
for case o of the two coupled oscillators. The finite, nonzero singular values are consistent with the surjectivity 
assumption in the analysis of Section [III 



state. Moreover, the Hessian at the optimal critical point has an infinite-dimensional nullspace, indicating 
an inherent degree of robustness to control field noise. The classical mechanical phase space state-to-state 
control results found here are clear analogues with the goal of maximizing quantum mechanical pure state 
transition probabilities. Interestingly, this analogy lies in correlating the number of classical particles with 
the number of quantum mechanical states. Further study is warranted to seek additional classical-quantum 
control analogies and relationships. 

Our findings extend those of [33J from open classical systems to closed, single trajectory systems. The 
analysis in [33J is kinematic, while the present paper started with the dynamical equations. The resulting 
analysis in Sections II A and II B rested on the nature of the objective functional and the regularity of 
the controls considered. These results have important fundamental implications for the control of classical 
mechanical systems with Hamiltonian dynamics, and in particular provide a firm foundation for controlled 
molecular dynamics of possibly complex polyatomic systems [551 ■ 



Numerical simulations confirmed the surjectivity of Sz{T)/Se{t), and no evidence was found for traps on the 
landscape in these studies. In practical molecular control applications, likely only a portion z' (t) of the state 
vector z{t) would be specified for control, implying that there would be a reduced demand on surjectivity for 
only 6z' {T)/6e{t). However, so-called singular controls, at which this set of derivatives is not surjective, also 
exist for both classical and quantum systems [351 112] • Such critical points are not yet well understood in the 
context of either classical or quantum optimal control, calling for additional analysis. Moreover, this paper 
examines only the control of single trajectory systems. Therefore, a natural next step to consider is that of 
controlling an ensemble of trajectories, with the state described by a probability distribution in the phase 
space in analogy with ensemble-based control of quantum systems [59j . Collectively, these studies may form 
a basis for considering control behavior that transcends the classical and quantum regimes. 
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APPENDIX: THE FORCED HARMONIC OSCILLATOR 



To illustrate the surjectivity assumption in Section II A consider the special case of a linear forced harmonic 



oscillator with potential V{q) 
see that 



(l/2)fcg^ and dipole moment D{q) ~ aq. Then, from Eqs. ^ and we 



A{t) = A 







B(t) = B 



which are both independent of time t. As a result, the matrix M{t) satisfying Eq. (12) can be simply 
expressed as M [t) = exp At, leading to the expression 



6e{t) 



= [cxpA{T -t)]B. 



We may explicitly determine M{t) by noting that 



A 



2Vfcr) 
1 
2 







^km 1 

V fcm 



For ease of notation, we define the harmonic oscillator frequency as w = \/{k/r 
A{T — t) and utilizing B{t) above yields 



Then exponentiating 



6e{t) 



exp [iuj{T - t)] + exp [iuj{t - T)] 
■.sm[uj{T -t)] 



/ km 

cos [cj(r - 1)] 

The two rows of Sz(T) / Se(t) are linear independent functions of time, thereby showing that the surjectivity 
assumption holds for this system. 

We further consider the special case of an anti-harmonic oscillator with potential V{q) = —{\/2)kq^ 
(k > 0) and dipole moment D(q) = aq. This system exhibits unstable behavior with z growing exponentially 
with time; nevertheless, it is a model for consideration of surjectivity analysis. As with molecular systems 
where dissociation occurs, the dynamics would only be followed for a finite time T out to when the particles 
are separated beyond an ability to impose further control. 

Following the forced harmonic oscillator above, we may derive 



Seit) 



7^ sinh \ujiT - t)] 

/km '- ^ '-^ 

cosh [lu{T - t)] 



The two rows of 62,{T) / Se{t) are linear independent functions of time, thereby showing that the surjectivity 
assumption holds for this extreme system. 



[1] C. Brif, R. Chakrabarti, and H. Rabitz, "Control of quantum phenomena: Past, present and future," New 
Journal of Physics, vol. 12, p. 075008, 2010. 

[2] H. Rabitz, R. de Vivie-Riedle, M. Motzkus, and K. Kompa, "Whither the future of controlling quantum phe- 
nomena?," Science, vol. 288, no. 5467, pp. 824-828, 2000. 

[3] W. S. Warren, H. Rabitz, and M. Dahleh, "Coherent control of quantum dynamics: The dream is alive," Science, 



18 



vol. 259, no. 5101, pp. 1581-1589, 1993. 

E. F. de Lima and M. A. M. do Aguiar, "Quantum-classical correspondence in the phase control of multiphoton 
dissociation by two-color laser pulses," Physical Review A, vol. 77, no. 3, p. 33406, 2008. 

R. B. Gerber, M. V. Korolkov, J. Manz, M. Y. Niv, and B. Schmidt, "A reflection principle for the control 
of molecular photodissociation in solids: Model simulation for F2 in Ar," Chemical Physics Letters, vol. 327, 
no. 1-2, pp. 76-84, 2000. 

S. Krcmpl, T. Eisenhammer, A. Hiibler, G. Mayer-Kress, and P. W. Milonni, "Optimal stimulation of a con- 
servative nonlinear oscillator: Classical and quantum-mechanical calculations," Physical Review Letters, vol. 69, 
pp. 430-433, Jul 1992. 

W. H. Miller, "A classical/semiclassical theory for the interaction of infrared radiation with molecular systems," 
The Journal of Chemical Physics, vol. 69, p. 2188, 1978. 

C. D. Schwieters and H. Rabitz, "Optimal control of nonlinear classical systems with application to unimolecular 

dissociation reactions and chaotic potentials," Physical Review A, vol. 44, no. 8, pp. 5224-5238, 1991. 

C. D. Schwieters and H. Rabitz, "Optimal control of classical systems with explicit quantum-classical-diffcrence 

reduction," Physical Review A, vol. 48, no. 4, pp. 2549-2557, 1993. 

R. B. Walker and R. K. Preston, "Quantum versus classical dynamics in the treatment of multiple photon 

excitation of the anharmonic oscillator," The Journal of Chemical Physics, vol. 67, pp. 2017-2028, 1977. 

V. S. Batista and P. Brumer, "Coherent control in the presence of intrinsic docoherence: Proton transfer in large 

molecular systems," Physical Review Letters, vol. 89, no. 14, p. 143201, 2002. 

J. Botina, H. Rabitz, and N. Rahman, "A new approach to molecular classical optimal control: Application to 
the reaction HON? HC-|- N," The Journal of Chemical Physics, vol. 102, pp. 226-236, 1995. 
J. Botina, H. Rabitz, and N. Rahman, "Optimal control of chaotic Hamiltonian dynamics," Physical Review A, 
vol. 51, no. 2, pp. 923-933, 1995. 

M. Demiralp and H. Rabitz, "Optimal control of cleissical molecular dynamics: A perturbation formulation and 

the existence of multiple solutions," Journal of Mathematical Chemistry, vol. 16, no. 1, pp. 185-209, 1994. 

A. Efimov, A. Fradkov, and A. Krivtsov, "Feedback design of control algorithms for dissociation of diatomic 

molecules," in Proceedings of the European Control Conference Cambridge, UK, pp. 1—4, 2003. 

V. Engel, C. Meier, and D. Tannor, "Applications to energy and particle transfer processes in molecules," in 

Advances in Chemical Physics (S. Rice, ed.), vol. 141 of Advances in Chemical Physics, John Wiley & Sons, 

2009. 

J. L. Krause, R. M. Whitnell, K. R. Wilson, and Y. J. Yan, ""Classical" quantum control with application to 

solution reaction dynamics," in AIP Conference Proceedings, vol. 298, pp. 3-15, Citeseer, 1993. 

M. H. Lissak, J. D. Sensabaugh, C. D. Schwieters, J. G. B. Bcumee, and H. Rabitz, "Optimal control of classical 

anharmonic molecules represented with picccwise harmonic potential surfaces: Analytic solution and selective 

dissociation of triatomic systems," The Journal of Chemical Physics, vol. 174, no. 1, pp. 1-24, 1993. 

Y. Nishiyama, T. Kato, Y. Ohtsuki, and Y. F^jimura, "Optimal laser control of ultrafast photodissociation of I 

in water: Mixed quantum/classical molecular dynamics simulation," The Journal of Chemical Physics, vol. 121, 

pp. 2685-2693, 2004. 

H. Umeda and Y. Fujimura, "Quantum control of chemical reaction dynamics in a classical way," The Journal 
of Chemical Physics, vol. 113, pp. 3510-3518, 2000. 

R. D. Taylor and P. Brumer, "Dynamical instability and external perturbations: Bimolecular collisions in laser 
fields," The Journal of Chemical Physics, vol. 77, p. 854, 1982. 

Y. J. Yan, R. E. Gillilan, R. M. Whitnell, K. R. Wilson, and S. Mukamel, "Optical control of molecular dynamics: 

Liouville-space theory," The Journal of Physical Chemistry, vol. 97, no. 10, pp. 2320- 2333, 1993. 

C. Cai, Z. Xu, and W. Xu, "Converting chaos into periodic motion by state feedback control* 1," Automatica, 

vol. 38, no. 11, pp. 1927-1933, 2002. 

Y. Lai, M. Ding, and C. Grebogi, "Controlling hamiltonian chaos," Physical Review E, vol. 47, no. 1, p. 86, 1993. 
E. Ott, C. Grebogi, and J. Yorke, "Controlling chaos," Physical review letters, vol. 64, no. 11, pp. 1196-1199, 
1990. 

Z. Wu, Z. Zhu, and C. Zhang, "Controlling hamiltonian chaos by medium perturbation in periodically driven 
systems," Physical Review E, vol. 57, no. 1, p. 366, 1998. 

T.-S. Ho and H. Rabitz, "Why do effective quantum controls appear easy to find?," Journal of Photochemistry 
and Photobiology A: Chemistry, vol. 180, no. 3, pp. 226-240, 2006. 

H. Rabitz, T.-S. Ho, M. Hsieh, R. Kosut, and M. Demiralp, "Topology of optimally controlled quantum mechan- 
ical transition probability landscapes," Physical Review A, vol. 74, no. 1, p. 12721, 2006. 

H. Rabitz, M. Hsieh, and C. Rosenthal, "Quantum optimally controlled transition landscapes," Science, vol. 303, 
no. 5666, pp. 1998-2001, 2004. 

M. Hsieh, R. Wu, and H. Rabitz, "Topology of the quantum control landscape for observables," The Journal of 
Chemical Physics, vol. 130, p. 104109, 2009. 



19 



R. Wu, H. Rabitz, and M. Hsieh, "Characterization of the critical submanifolds in quantum ensemble control 
landscapes," Journal of Physics A: Mathematical and Theoretical, vol. 41, p. 015006, 2008. 

R. Wu, R. Chakrabarti, and H. Rabitz, "Critical landscape topology for optimization on the symplectic group," 
Journal of Optimization Theory and Applications, vol. 145, no. 2, pp. 387-406, 2010. 

A. Pechen and H. Rabitz, "Unified analysis of terminal-time control in classical and quantum systems," Euro- 
physics Letters, vol. 91, p. 60005, 2010. 

R. V. Mendes and V. I. Manko, "Quantum control and the Strocchi map," Physical Review A, vol. 67, no. 5, 
p. 053404, 2003. 

A. Heslot, "Quantum mechanics as a classical theory," Physical Review D, vol. 31, no. 6, pp. 1341-1348, 1985. 
Y. M. Shirokov, "Quantum and classical mechanics in the phase space representation," Soviet Journal of Nuclear 
Physics, vol. 10, no. 1, pp. 1-18, 1979. 

E. D. Sontag, Mathematical Control Theory: Deterministic Finite Dimensional Systems. Texts in applied math- 
ematics. Springer, 1998. 

R. Abrines and I. Percival, "Classical theory of charge transfer and ionization of hydrogen atoms by protons," 
Proceedings of the Physical Society, vol. 88, p. 861, 1966. 

B. Bonnard and M. Chyba, Singular Trajectories and Their Role in Control Theory. Springer Verlag, 2003. 
V. Jurdjevic, Geometric Control Theory. New York: Cambridge UP, 1997. 

Sufficient conditions are known under which the reachable set has a non-empty interior. However, it is difhcult 
to derive general necessary and sufficient conditions [40] . 

R. Wu, J. Dominy, T.-S. Ifo, and H. Rabitz, "Singularities of quantum control landscapes," Arxiv preprint, 2009. 
arXiv:0907.2354. 

R. W. Brockett, Finite Dimensional Linear Systems. Series in decision and control, Wiley, 1970. 

A. J. Dragt, "The symplectic group and classical mechanics," Annals of the New York Academy of Sciences, 

vol. 1045, no. 1, pp. 291-307, 2005. 

K. Hoffman and R. A. Kunze, Linear Algebra. Prentice-Hall mathematics series, Prentice-Hall, 1971. 

T. Ho, J. Dominy, and H. Rabitz, "Landscape of unitary transformations in controlled quantum dynamics," 

Physical Review A, vol. 79, no. 1, p. 013422, 2009. 

A. Rothman, T.-S. Ho, and H. Rabitz, "Exploring the level sets of quantum control landscapes," Physical Review 
A, vol. 73, no. 5, p. 53401, 2006. 

V. Beltrani, J. Dominy, T.-S. Ho, and H. Rabitz, "Exploring the top and bottom of the quantum control 
landscape," Arxiv preprint, 2011. arXiv:1103.5390. 

D. Donnelly and E. Rogers, "Symplectic integrators: An introduction," American Journal of Physics, vol. 73, 
pp. 938-945, 2005. 

A. Rothman, T.-S. Ho, and H. Rabitz, "Observable-preserving control of quantum dynamics over a family of 
related systems," Physical Review A, vol. 72, no. 2, p. 23416, 2005. 

J. Candy and W. Rozmus, "A symplectic integration algorithm for separable Hamiltonian functions," Journal 
of Computational Physics, vol. 92, no. 1, pp. 230-256, 1991. 

R. I. McLachlan and P. Atela, "The accuracy of symplectic integrators," Nonlinearity, vol. 5, pp. 541-562, 1992. 
A. Rothman, T.-S. Ho, and H. Rabitz, "Quantum observable homotopy tracking control," The Journal of Chem- 
ical Physics, vol. 123, p. 134104, 2005. 

A. Guldberg and G. D. Billing, "Laser-induced dissociation of hydrogen fluoride," Chemical Physics Letters, 
vol. 186, no. 2-3, pp. 229-237, 1991. 

J. R. Stine and D. W. Noid, "Classical treatment of the dissociation of hydrogen fluoride with one and two 
infrared lasers," Optics Communications, vol. 31, no. 2, pp. 161-164, 1979. 

M. L. Sage and J. A. Williams HI, "Energetics, wave functions, and spectroscopy of coupled anharmonic oscil- 
lators," The Journal of Chemical Physics, vol. 78, p. 1348, 1983. 

M. D. Radicioni, C. G. Diaz, and F. M. Fernandez, "Application of perturbation theory to coupled morse 
oscillators," Journal of Molecular Structure: THEOCHEM, vol. 488, no. 1-3, pp. 37-49, 1999. 
D. C. Rapaport, The Art of Molecular Dynamics Simulation. Cambridge UP, 2004. 

H. Rabitz, M. Hsieh, and C. Rosenthal, "Optimal control landscapes for quantum observables," The Journal of 
Chemical Physics, vol. 124, p. 204107, 2006. 



20 



